/***************************************************************************\
 * Copyright (C) by Keio University
 * BilateralFilter.cu created in 01 2012.
 * Mail : fdesorbi@hvrl.ics.keio.ac.jp
 *
 * BilateralFilter.cu is part of the HVRL Engine Library.
 *
 * The HVRL Engine Library is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * The HVRL Engine Library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 \***************************************************************************/

#ifndef BILATERAL_FILTER_CU_
#define BILATERAL_FILTER_CU_

#ifdef WIN32
#define NOMINMAX
#define M_PI 3.141596
#include <windows.h>
#endif

#include <hvrl/gpgpu/cuda/Prototypes.hpp>
#include <hvrl/gpgpu/cuda/Tools.cuh>

namespace hvrl {

namespace cuda {

namespace BilateralFilter {

#define MAX_KERNEL_RADIUS 8
#define MAX_KERNEL_LENGTH (2 * MAX_KERNEL_RADIUS + 1)
#define MAX_KERNEL_SIZE MAX_KERNEL_LENGTH*MAX_KERNEL_LENGTH

__device__ __constant__ float d_spatial_kernel[MAX_KERNEL_SIZE]; /**< constant size array for holding the spatial kernel located in the constant memory of the host*/
__device__ __constant__ int d_radius; /**< Radius of the kernel located in the constant memory of the host*/
__device__ __constant__ float d_coefficient_range; /**< Precomputed coefficient for the range filter located in the constant memory of the host*/
__device__ __constant__ float d_divider_range; /**< Precomputed coefficient for the range filter divider located in the constant memory of the host*/

/**
 * \struct Kernel
 * \brief Holds the parameters related to the bilateral filter kernel plus some other pre-processed values
 * \author Francois de Sorbier
 */
struct Kernel {
	/**
	 * \name internal parameters of the bilateral filter
	 */
	/*\{*/
	int radius; /**< Radius of the kernel */
	/*\}*/

	/**
	 * \name pre-computed values
	 */
	/*\{*/
	float coefficient_range; /**< Hold the pre-computed the coefficient 1 / (sqrt(2*PI)*sigma) */
	float divider_range; /**< Hold the pre-computed the coefficient - 1 / (2*sigma^2) */
	float spatial_kernel[MAX_KERNEL_SIZE]; /**< Hold the pre-computed spatial kernel on the device */
	/*\}*/
};

/**
 * \brief Allocate and fill the structure for the kernel
 *
 * The spatial kernel is computed as \f$G_{\sigma_s}(x,y) = \frac{1}{2\pi\sigma_s^2}\exp({{-}{x^2+y^2}\over{2\sigma_s^2}})\f$
 *
 * \param[in] sigma_s the spatial gaussian coefficient
 * \param[in] sigma_r the range gaussian coefficient
 * \param[in] radius of the kernel. The maximum size is defined by MAX_KERNEL_RADIUS for better efficiency
 * \return the bilateral kernel allocated on the device
 * \author Francois de Sorbier
 */
Kernel*
getKernel(const float& sigma_s, const float& sigma_r, const int& radius) {

	// We need a copy of the structure on the host
	Kernel* output = new Kernel();

	float coefficient_space = 1.0f / (2.0f * M_PI * powf(sigma_s, 2));
	output->coefficient_range = 1.0f / (sqrtf(2.0f * M_PI) * sigma_r);
	float divider_space = -1.0f / (2.0f * powf(sigma_s, 2));
	output->divider_range = -0.5f / (powf(sigma_r, 2));
	if (radius > MAX_KERNEL_RADIUS)
		output->radius = MAX_KERNEL_RADIUS;
	else
		output->radius = radius;

	int width = output->radius * 2 + 1;

	// Precompute the 2D spatial kernel
	for (int i = -output->radius; i <= output->radius; ++i) {
		for (int j = -output->radius; j <= output->radius; ++j) {
			output->spatial_kernel[(i + output->radius) * width
					+ (j + output->radius)] = coefficient_space
					* exp(
							(pow(float(i), 2) + pow(float(j), 2))
									* divider_space);
		}
	}

	return output;
}

/**
 * \brief Create and precompute parts of the BF kernel on the device
 *
 * The spatial kernel is precomputed as \f$G_{\sigma_s}(x,y) = \frac{1}{2\pi\sigma_s^2}\exp({{-}{x^2+y^2}\over{2\sigma_s^2}})\f$
 *
 * \param[in] sigma_s the spatial gaussian coefficient
 * \param[in] sigma_r the range gaussian coefficient
 * \param[in] radius of the kernel. The maximum size is defined by MAX_KERNEL_RADIUS for better efficiency
 * \sa getKernel setKernel
 * \author Francois de Sorbier
 */
void setKernel(const float& sigma_s, const float& sigma_r, const int& radius) {

	float coefficient_space = 1.0f / (2.0f * M_PI * powf(sigma_s, 2));
	float coefficient_range = 1.0f / (sqrtf(2.0f * M_PI) * sigma_r);
	float divider_space = -1.0f / (2.0f * powf(sigma_s, 2));
	float divider_range = -0.5f / (powf(sigma_r, 2));
	int tmp_radius = MAX_KERNEL_RADIUS;
	if (radius < MAX_KERNEL_RADIUS && radius > 0)
		tmp_radius = radius;

	checkStatus(cudaMemcpyToSymbol(d_radius, &tmp_radius, sizeof(int)));
	checkStatus(
			cudaMemcpyToSymbol(d_coefficient_range, &coefficient_range,
					sizeof(float)));
	checkStatus(
			cudaMemcpyToSymbol(d_divider_range, &divider_range, sizeof(float)));

	int width = tmp_radius * 2 + 1;
	float kernel[MAX_KERNEL_SIZE];

	// Precompute the 2D spatial kernel
	for (int i = -tmp_radius; i <= tmp_radius; ++i) {
		for (int j = -tmp_radius; j <= tmp_radius; ++j) {
			kernel[(i + tmp_radius) * width + (j + tmp_radius)] =
					coefficient_space
							* exp(
									(pow(float(i), 2) + pow(float(j), 2))
											* divider_space);
		}
	}

	checkStatus(
			cudaMemcpyToSymbol(d_spatial_kernel, kernel,
					MAX_KERNEL_SIZE * sizeof(float)));
}

/**
 * \brief Create the BF kernel on the device based on a precomputed Kernel from the host
 *
 * \param[in] kernel the precomputed kernel
 * \sa getKernel setKernel
 * \author Francois de Sorbier
 */
void setKernel(Kernel* kernel) {
	checkStatus(cudaMemcpyToSymbol(d_radius, &(kernel->radius), sizeof(int)));
	checkStatus(
			cudaMemcpyToSymbol(d_coefficient_range,
					&(kernel->coefficient_range), sizeof(float)));
	checkStatus(
			cudaMemcpyToSymbol(d_divider_range, &(kernel->divider_range),
					sizeof(float)));
	checkStatus(
			cudaMemcpyToSymbol(d_spatial_kernel, kernel->spatial_kernel,
					MAX_KERNEL_SIZE * sizeof(float)));
}

/**
 * \brief The Kernel computes the result of a bilateral filtering (approximated) for a grayscale image
 *
 * Bilateral filter :
 *
 * \f$I_p = \frac{1}{W_p} \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)I_q}\f$
 *
 * with \f$W_p = \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)}\f$
 *
 * \f$G_{\sigma}(x) = \exp({{-}{x^2}\over{2\sigma^2}})\f$
 *
 * \f$G_{\sigma_s}\f$ is named the space weight (distance from center) and \f$G_{\sigma_r}\f$ is the range weight (difference of intensity)
 *
 * In our case, \f$G_{\sigma_s}(\|p-q\|)\f$ is approximated to \f$G_{\sigma_s}(\|p_x-q_x\|)G_{\sigma_s}(\|p_y-q_y\|)\f$
 *
 * \param output contains the result of the bilateral filtering
 * \param input the input grayscale image
 * \param e_d Gaussian range parameter \f$ 0.1 \le \sigma_r < \infty \f$
 * \param w the width of the input/output
 * \param h the height of the input/output
 * \param r radius of the filter
 * \author Francois de Sorbier
 */
__global__ void d_bilateral_filter_gray(unsigned char *output,
		unsigned char *input, const int w, const int h) {
	int x = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
	int y = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;

	if (x < w && x >= 0 && y < h && y >= 0) {
		float sum = 0.0f;
		float t = 0.0f;

		// Need to convert to normalized float
		float center = float(input[y * w + x] * 0.003921569f);

		for (int i = -d_radius; i <= d_radius; i++) {
			if (y + i >= 0 && y + i < h) {
				for (int j = -d_radius; j <= d_radius; j++) {
					// Check we are in the limits and in the circle
					if (x + j >= 0 && x + j < w) {

						float curPix = float(
								input[(y + i) * w + x + j] * 0.003921569f);
						float eulen = center - curPix;
						float g = d_spatial_kernel[(i + d_radius)
								* (d_radius * 2 + 1) + j + d_radius]
								* d_coefficient_range
								* __expf(pow(eulen, 2) * d_divider_range);
						t += g * curPix;
						sum += g;
					}
				}
			}

		}
		output[y * w + x] = (unsigned char) (t / sum * 255.0f);
	}
}

/**
 * \brief The Kernel computes the result of a bilateral filtering (approximated) for a grayscale image
 *
 * Bilateral filter :
 *
 * \f$I_p = \frac{1}{W_p} \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)I_q}\f$
 *
 * with \f$W_p = \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)}\f$
 *
 * \f$G_{\sigma}(x) = \exp({{-}{x^2}\over{2\sigma^2}})\f$
 *
 * \f$G_{\sigma_s}\f$ is named the space weight (distance from center) and \f$G_{\sigma_r}\f$ is the range weight (difference of intensity)
 *
 * In our case, \f$G_{\sigma_s}(\|p-q\|)\f$ is approximated to \f$G_{\sigma_s}(\|p_x-q_x\|)G_{\sigma_s}(\|p_y-q_y\|)\f$
 *
 * \param output contains the result of the bilateral filtering
 * \param input the input grayscale image
 * \param e_d Gaussian range parameter \f$ 0.1 \le \sigma_r < \infty \f$
 * \param w the width of the input/output
 * \param h the height of the input/output
 * \param r radius of the filter
 * \author Francois de Sorbier
 */
__global__ void d_bilateral_filter_gray(short unsigned int *output,
		short unsigned int *input, const int w, const int h) {
	int x = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
	int y = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;

	if (x < w && x >= 0 && y < h && y >= 0) {
		float sum = 0.0f;
		float t = 0.0f;

		// Need to convert to normalized float
		float center = float(input[y * w + x]);

		for (int i = -d_radius; i <= d_radius; i++) {
			if (y + i >= 0 && y + i < h) {
				for (int j = -d_radius; j <= d_radius; j++) {
					// Check we are in the limits and in the circle
					if (x + j >= 0 && x + j < w) {

						float curPix = float(input[(y + i) * w + x + j]);
						float eulen = center - curPix;
						float g = d_spatial_kernel[(i + d_radius)
								* (d_radius * 2 + 1) + j + d_radius]
								* d_coefficient_range
								* __expf(pow(eulen, 2) * d_divider_range);
						t += g * curPix;
						sum += g;
					}
				}
			}

		}
		output[y * w + x] = (short) (t / sum);
	}
}

/**
 * \brief The Kernel computes the result of a bilateral filtering (approximated) for a grayscale image
 *
 * Bilateral filter :
 *
 * \f$I_p = \frac{1}{W_p} \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)I_q}\f$
 *
 * with \f$W_p = \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)}\f$
 *
 * \f$G_{\sigma}(x) = \exp({{-}{x^2}\over{2\sigma^2}})\f$
 *
 * \f$G_{\sigma_s}\f$ is named the space weight (distance from center) and \f$G_{\sigma_r}\f$ is the range weight (difference of intensity)
 *
 * In our case, \f$G_{\sigma_s}(\|p-q\|)\f$ is approximated to \f$G_{\sigma_s}(\|p_x-q_x\|)G_{\sigma_s}(\|p_y-q_y\|)\f$
 *
 * \param output contains the result of the bilateral filtering
 * \param input the input grayscale image
 * \param e_d Gaussian range parameter \f$ 0.1 \le \sigma_r < \infty \f$
 * \param w the width of the input/output
 * \param h the height of the input/output
 * \param r radius of the filter
 * \author Francois de Sorbier
 */
__global__ void d_bilateral_filter_gray(float *output, float *input,
		const int w, const int h) {
	int x = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
	int y = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;

	if (x < w && x >= 0 && y < h && y >= 0) {
		float sum = 0.0f;
		float t = 0.0f;

		// Need to convert to normalized float
		float center = input[y * w + x];

		for (int i = -d_radius; i <= d_radius; i++) {
			if (y + i >= 0 && y + i < h) {
				for (int j = -d_radius; j <= d_radius; j++) {
					// Check we are in the limits and in the circle
					if (x + j >= 0 && x + j < w) {

						float curPix = input[(y + i) * w + x + j];
						float eulen = center - curPix;
						float g = d_spatial_kernel[(i + d_radius)
								* (d_radius * 2 + 1) + j + d_radius]
								* d_coefficient_range
								* __expf(pow(eulen, 2) * d_divider_range);
						t += g * curPix;
						sum += g;
					}
				}
			}

		}
		output[y * w + x] = t / sum;
	}
}

/**
 * \brief The Kernel computes the result of a bilateral filtering (approximated) for a color image
 *
 * Bilateral filter :
 *
 * \f$I_p = \frac{1}{W_p} \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)I_q}\f$
 *
 * with \f$W_p = \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)}\f$
 *
 * \f$G_{\sigma}(x) = \exp({{-}{x^2}\over{2\sigma^2}})\f$
 *
 * \f$G_{\sigma_s}\f$ is named the space weight (distance from center) and \f$G_{\sigma_r}\f$ is the range weight (difference of intensity)
 *
 * In our case, \f$G_{\sigma_s}(\|p-q\|)\f$ is approximated to \f$G_{\sigma_s}(\|p_x-q_x\|)G_{\sigma_s}(\|p_y-q_y\|)\f$
 *
 * \param output contains the result of the bilateral filtering
 * \param input the input color image
 * \param e_d Gaussian range parameter \f$ 0.1 \le \sigma_r < \infty \f$
 * \param w the width of the input/output
 * \param h the height of the input/output
 * \param r radius of the filter
 * \todo the computed could be speed-up by converting the color image to greyscale before hand or passing it as an argument
 * \author Francois de Sorbier
 */
__global__ void d_bilateral_filter_color(unsigned char *output,
		unsigned char *input, const int w, const int h) {
	int x = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
	int y = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;

	float3 curPix, center, t;
	float center_gray, curPix_gray, sum;

	if (x < w && x >= 0 && y < h && y >= 0) {
		sum = 0.0f;
		t.x = t.y = t.z = 0.0f;

		// Need to convert to normalized float
		center.x = float(input[(y * w + x) * 3] * 0.003921569f);
		center.y = float(input[(y * w + x) * 3 + 1] * 0.003921569f);
		center.z = float(input[(y * w + x) * 3 + 2] * 0.003921569f);
		center_gray = center.x * 0.3f + center.y * 0.59f + center.z * 0.11f;

		for (int i = -d_radius; i <= d_radius; i++) {
			if (y + i >= 0 && y + i < h) {
				for (int j = -d_radius; j <= d_radius; j++) {
					// Check we are in the limits and in the circle
					if (x + j >= 0 && x + j < w) {

						curPix.x = float(
								input[((y + i) * w + x + j) * 3]
										* 0.003921569f);
						curPix.y = float(
								input[((y + i) * w + x + j) * 3 + 1]
										* 0.003921569f);
						curPix.z = float(
								input[((y + i) * w + x + j) * 3 + 2]
										* 0.003921569f);
						curPix_gray = curPix.x * 0.3f + curPix.y * 0.59f
								+ curPix.z * 0.11f;
						float eulen = center_gray - curPix_gray;
						float g = d_spatial_kernel[(i + d_radius)
								* (d_radius * 2 + 1) + j + d_radius]
								* d_coefficient_range
								* __expf(pow(eulen, 2) * d_divider_range);
						t.x += g * curPix.x;
						t.y += g * curPix.y;
						t.z += g * curPix.z;
						sum += g;
					}
				}
			}

		}
		sum = 255.0f / sum;
		output[(y * w + x) * 3] = (unsigned char) (t.x * sum);
		output[(y * w + x) * 3 + 1] = (unsigned char) (t.y * sum);
		output[(y * w + x) * 3 + 2] = (unsigned char) (t.z * sum);
	}
}

/**
 * \brief The Kernel computes the result of a bilateral filtering (approximated) for a color image
 *
 * Bilateral filter :
 *
 * \f$I_p = \frac{1}{W_p} \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)I_q}\f$
 *
 * with \f$W_p = \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)}\f$
 *
 * \f$G_{\sigma}(x) = \exp({{-}{x^2}\over{2\sigma^2}})\f$
 *
 * \f$G_{\sigma_s}\f$ is named the space weight (distance from center) and \f$G_{\sigma_r}\f$ is the range weight (difference of intensity)
 *
 * In our case, \f$G_{\sigma_s}(\|p-q\|)\f$ is approximated to \f$G_{\sigma_s}(\|p_x-q_x\|)G_{\sigma_s}(\|p_y-q_y\|)\f$
 *
 * \param output contains the result of the bilateral filtering
 * \param input the input color image
 * \param e_d Gaussian range parameter \f$ 0.1 \le \sigma_r < \infty \f$
 * \param w the width of the input/output
 * \param h the height of the input/output
 * \param r radius of the filter
 * \todo the computed could be speed-up by converting the color image to greyscale before hand or passing it as an argument
 * \author Francois de Sorbier
 */
__global__ void d_bilateral_filter_color(short unsigned int *output,
		short unsigned int *input, const int w, const int h) {
	int x = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
	int y = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;

	float3 curPix, center, t;
	float center_gray, curPix_gray, sum;

	if (x < w && x >= 0 && y < h && y >= 0) {
		sum = 0.0f;
		t.x = t.y = t.z = 0.0f;

		// Need to convert to normalized float
		center.x = float(input[(y * w + x) * 3]);
		center.y = float(input[(y * w + x) * 3 + 1]);
		center.z = float(input[(y * w + x) * 3 + 2]);
		center_gray = center.x * 0.3f + center.y * 0.59f + center.z * 0.11f;

		for (int i = -d_radius; i <= d_radius; i++) {
			if (y + i >= 0 && y + i < h) {
				for (int j = -d_radius; j <= d_radius; j++) {
					// Check we are in the limits and in the circle
					if (x + j >= 0 && x + j < w) {

						curPix.x = float(input[((y + i) * w + x + j) * 3]);
						curPix.y = float(input[((y + i) * w + x + j) * 3 + 1]);
						curPix.z = float(input[((y + i) * w + x + j) * 3 + 2]);
						curPix_gray = curPix.x * 0.3f + curPix.y * 0.59f
								+ curPix.z * 0.11f;
						float eulen = center_gray - curPix_gray;
						float g = d_spatial_kernel[(i + d_radius)
								* (d_radius * 2 + 1) + j + d_radius]
								* d_coefficient_range
								* __expf(pow(eulen, 2) * d_divider_range);
						t.x += g * curPix.x;
						t.y += g * curPix.y;
						t.z += g * curPix.z;
						sum += g;
					}
				}
			}

		}
		sum = 1.0f / sum;
		output[(y * w + x) * 3] = (short) (t.x * sum);
		output[(y * w + x) * 3 + 1] = (short) (t.y * sum);
		output[(y * w + x) * 3 + 2] = (short) (t.z * sum);
	}
}

/**
 * \brief Launch the kernel that will apply the bilateral filter for a grayscale image
 * \param h_input the input grayscale image
 * \param h_output the result of the bilateral filtering
 * \param h_width the width of the input image
 * \param h_height the height of the output image
 * \author Francois de Sorbier
 */
void computeC1(unsigned char *d_input, unsigned char *d_output,
		const int& d_width, const int& d_height) {
	// Define the size of the grid and the block
	dim3 blockSize(32, 32);
	dim3 gridSize((d_width + blockSize.x - 1) / blockSize.x,
			(d_height + blockSize.y - 1) / blockSize.y);

d_bilateral_filter_gray <<< gridSize, blockSize>>>(d_output, d_input, d_width, d_height);

}

/**
 * \brief Launch the kernel that will apply the bilateral filter for a color image
 * \param h_input the input color image
 * \param h_output the result of the bilateral filtering
 * \param h_width the width of the input image
 * \param h_height the height of the output image
 * \author Francois de Sorbier
 */
void computeC3(unsigned char *d_input, unsigned char *d_output,
		const int& d_width, const int& d_height) {
	// Define the size of the grid and the block
	dim3 blockSize(32, 32);
	dim3 gridSize((d_width + blockSize.x - 1) / blockSize.x,
			(d_height + blockSize.y - 1) / blockSize.y);

d_bilateral_filter_color <<< gridSize, blockSize>>>(d_output, d_input, d_width, d_height);

}

/**
 * \brief Launch the kernel that will apply the bilateral filter for a grayscale image
 * \param h_input the input grayscale image
 * \param h_output the result of the bilateral filtering
 * \param h_width the width of the input image
 * \param h_height the height of the output image
 * \author Francois de Sorbier
 */
void computeC1(short unsigned int *d_input, short unsigned int *d_output,
		const int& d_width, const int& d_height) {
	// Define the size of the grid and the block
	dim3 blockSize(32, 32);
	dim3 gridSize((d_width + blockSize.x - 1) / blockSize.x,
			(d_height + blockSize.y - 1) / blockSize.y);

d_bilateral_filter_gray <<< gridSize, blockSize>>>(d_output, d_input, d_width, d_height);

}

/**
 * \brief Launch the kernel that will apply the bilateral filter for a grayscale image
 * \param h_input the input grayscale image
 * \param h_output the result of the bilateral filtering
 * \param h_width the width of the input image
 * \param h_height the height of the output image
 * \author Francois de Sorbier
 */
void computeC1(float *d_input, float *d_output, const int& d_width,
		const int& d_height) {
	// Define the size of the grid and the block
	dim3 blockSize(32, 32);
	dim3 gridSize((d_width + blockSize.x - 1) / blockSize.x,
			(d_height + blockSize.y - 1) / blockSize.y);

d_bilateral_filter_gray <<< gridSize, blockSize>>>(d_output, d_input, d_width, d_height);

}

/**
 * \brief Launch the kernel that will apply the bilateral filter for a color image
 * \param h_input the input color image
 * \param h_output the result of the bilateral filtering
 * \param h_width the width of the input image
 * \param h_height the height of the output image
 * \author Francois de Sorbier
 */
void computeC3(short unsigned int *d_input, short unsigned int *d_output,
		const int& d_width, const int& d_height) {
	// Define the size of the grid and the block
	dim3 blockSize(32, 32);
	dim3 gridSize((d_width + blockSize.x - 1) / blockSize.x,
			(d_height + blockSize.y - 1) / blockSize.y);

d_bilateral_filter_color <<< gridSize, blockSize>>>(d_output, d_input, d_width, d_height);

}

}



}
}

#endif
